Thermodynamics of dissipative quantum systems by effective potential 
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Classical-like formulas are given in order to evaluate thermal averages of observables belonging 
to a quantum nonlinear system with dissipation described by the Caldeira-Leggett model [Phys. 
Rev. Lett. 46, 211 (1981); Ann. Phys. (N.Y.) 149, 374 (1983)]. The underlying scheme is the 
pure-quantum self-consistent harmonic approximation, which leads to expressions with a Boltzmann 
factor involving an effective potential and with a Gaussian average. The latter describes the effect of 
the fluctuations of purely quantum origin. As an illustration we calculate the coordinate probability 
distribution for a double-well potential in the presence of various degrees of Ohmic dissipation. 



The concept of effective potential, meant to reduce 
quantum statistical mechanics calculations to classical 
ones was firs introduced by Feynman Jl],|| . He introduced 
a variational principle for the path-integral expression of 
the partition function — the Feynman-. J ens en (FJ) in- 
equality — and used it with a "free particle" trial action. 

A significant improvement has been achieved by Gi- 
achetti and Tognetti and independently by Feyn- 
man and Kleinert || using a quadratic trial action with 
the same variational principle. For nonstandard sys- 
tems, where the FJ inequality is generally not valid, 
the pure-quantum self- consistent harmonic approxima- 
tion (PQSCHA) gives a way to construct an effective 
Hamiltonian, thus recovering the phase-space concept 
and the classical- like formulas for thermal averages ||[7| . 

Several successful applications in different branches of 
condensed matter physics pointed out the power of the 
approach 0. In the case of open systems, little work 
has been done for taking into account the quantum dis- 
sipation in the effective potential formalism; indeed, in 
a system-plus-reservoir model, only the expression of the 
effective potential for expressing the partition function 
as a configuration integral has been given, both for lin- 
ear and nonlinear coupling JlC| ] with environmental 
oscillators. 

By using the PQSCHA approach, that is equivalent 
to the variational method, we present here an accurate 
derivation of the density matrix of a nonlinear system 
interacting with a heat bath through the Caldeira-Leggett 
(CL) model jnj. This model considers the system of 
interest as linearly interacting with a bath of harmonic 
oscillators, whose coordinates can be integrated out from 
the path integral, leaving the CL Euclidean action 
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the environmental bath |9| . The density matrix elements 
in the coordinate representation are expressed by Feyn- 
man's path integral as 



p(q",q')={q"\p\ q ') 



V[q{u)\ e~ S ^ , (2) 



where the path integration is defined as a sum over all 
paths q(u), with u G [0, (3%], g(0) = q' , and q{f3K) = q" . 

As suggested by Feynman M, we can rearrange the 
path integral (Q) summing over classes of paths that share 
the same average point 



P(q",q') = J dqp(q",q';q) , 
p(q",q';q)= f V[q(u)]S(q - q[q(u)}) e 



-S[q[u)] 
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where q[q(u)] = (f3h)^ 1 J Q h duq(u) is the average-point 
functional. Furthermore, as only paths with a fixed av- 
erage point q appear into the path integral (|]) , from ac- 
tion (|l|) we define a trial action So^u)] by replacing 
V(q(uyj with a trial quadratic "potential" 



V (q; q) = w(q) + \muj 2 {q) (q - q) 2 



(5) 



The kernel k{u) = k(J3h — u) is a function that depends 
on temperature T = 1/(3 and on the spectral density of 



where the parameters w{q) and uj 2 (q) are now to be opti- 
mized, so that the trial reduced density po{q, q; q) at best 
approximates p(q,q;q), for each value of q. Note that 
Vo(q;q) is not related to a quantum observable, since it 
depends on q: so, the first part of the trial action is also 
nonlocal. 

The evaluation of the trial reduced density po(q" , q'',q), 
in spite of the fact that the action is quadratic, is rather 
tricky: on the one hand, the known general method of 
calculating the minimal action fails since its minimization 
gives rise to infinite order equations of motion, and, on 
the other hand, the method of Fourier expansion of the 
paths in terms of discrete Matsubara components con- 
flicts with the openness of the paths, since we are indeed 
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looking also for the off-diagonal part of the density ma- 
trix (q" ^ q'). However, the second way may still be 
followed by transforming as 



P(q",q';q) = 

lun | y j> T>[q(u)} s(g- q[q{u)} 



-S[q(u)} 



(6) 



the integral is over all closed paths {q{u) \ u S [0, (3h] } 
that satisfy the constraints q(e) = q' and q(f3% — e) = 
q" . T E is the integral over the open paths { qiu) \ u £ 
[— e, e] } (the range [(3K — e, (3K\ is periodically mapped 
onto [— e, 0]) with end points q(— e) = q" and q(e) = q'\ 
for small e it is dominated by the kinetic contribution 



T, = 
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Now the paths in Eq. ([|) can be Fourier expanded, 

oo 

q(u) = q + 2 2J (x n cos v n u + y n sin v n u) ; 



(7) 



(8) 



v n = 2Tm/((3h) are the Matsubara frequencies and the 
vq component is just q, as it appears from the inverse 
transformation. The measure of the path integral then 
becomes Hi 
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and the trial action takes then the form 
So[q{u)] = (3w(q) 

oc 

+f3m ^ [vl + u 2 (q) + VnliUnj] {x 2 n 
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Here we have made use of the relation that connects the 
kernel k(u) with the Laplace transform j(z) of the real- 
time memory damping function j(t) H 
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While it is trivial to manage the S function that fixes q, 
the end-point constraints are implemented by inserting 
the 5 functions 8(q(j3Ti—e) — q") and 6(q(e)—q'), and then 
using their Fourier representation. The calculation of 
po can then be carried forward by Gaussian quadratures 
giving 
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where £ = \{q' + q") — q and ( = q" — q' , 
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In the last expressions, the leading terms for small e are 
c E = 2a(q) + 0(e) , and, in a less straightforward way 
s e = 4he/m[l - ^X(q) + o(e)] , with 
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Pm ^ K + uj z (q) + v nl (v n ) 



v 2 n + ^ 2 (q) + Kh(Wn\) 
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(17) 



Eventually, the limit in e can be easily taken, and the 
result reads 



po{q",q';q) 
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Now, if Po(q",q') - obtained using the last result in 
Eq. ([?]) - is taken as an approximation for the exact den- 
sity matrix p(q" ,q'), we have a way to explicitly calculate 
any quantum thermal average by two Gaussian quadra- 
tures plus a single integral in q. A more convenient for- 
malism deals with phase space: let us briefly show it. 

The Weyl symbol for an observable O — 0(p, q) is the 
phase-space function 0(jp,q) that is defined in terms of 
the matrix elements in coordinate space as [fl3| 



0(p,«)= / d( e-^/ R (<z+iC|<%-K> 



(19) 



By means of a simple property of Weyl symbols the av- 
erage of O takes the form of a phase-space integral, 
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From Eqs. (|3|) and ( |18|) the Weyl symbol for the trial 
density operator turns out to be 



Po(p,<l) = ZkH 



2nh 2 p 



dq ■ 



p(q) 



e -S72a(g) e -V/2A(?) 



\/27ra(g) ^/2nX(q) 



(21) 



where £ = q — q. Therefore, the average of any observable 
0(p, q) can be expressed as a classical formula 



(6) 
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dq(O(p,q + 0} e"^® , (22) 
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where the double bracket is the Gaussian average op- 
erating over p and £, with moments ((£ 2 )) = a(q) and 
((p 2 )) = X(q); the effective potential is defined as 



V eS (q) = w(q) + (3 1 ]n f i(q) 



(23) 



In order to determine the functions w(q) and 0J 2 {q) 
we impose the PQSCHA condition, i.e., we require that 
the potential V(q) and the trial potential Vo(q;q), to- 
gether with their derivatives up to second order, have 
equal averages with respect to the reduced (diagonal) 
density po(q, q; q). The condition for the first derivatives 
is overcome by the definition of the average point, so that 
we are left with 

(V(q+0) = (TO+O) = + ~ ^(o) «(?) , 

(V"(q+0) = (V{'(q+Q) = mco 2 (q) . (24) 

The latter condition must be solved self-consistently with 
the definition of a{q) in terms of LD 2 (q), Eq. (|l6|). Now 
we have all the necessary ingredients to explicitly eval- 
uate the effective potential and all the thermal averages 
through the classical- like expression ( |22| ) . As in the case 
of no dissipation, it can be seen that for the most usual 
potentials the self-consistent solution for a(q) turns out 
to be always positive, even though u> 2 (q) can be negative 
HI). Indeed, what matters is that a, Eq. fll6]), is a de- 
creasing function of uj 2 , with a divergence to +oo [which 
in the dissipative case happens at uj 2 — > — v\ — vyyiyx) ]. 

In the case of Ohmic dissipation, corresponding to 
7(z) = 7 =const [i.e., the memory is Markovian, j(t) = 
7#(t — 0)], one can see that the (g-dependent) mean- 
square momentum X(q), Eq. (|l7j), is divergent in the 
Ohmic case. Correspondingly, also /«(<?), Eq. ( [l3| ) di- 
verges. The physical reason for this, basically related to 
the uncertainty principle, is well discussed in Refs. P,|lT||. 
Here we note that the coordinate mean-square fluctu- 
ation a(q) is still well defined, and that the effective 
potential can be made finite by subtraction of the in- 
finite but g-independent quantity j3~ lrxMi, with /zi = 
n^Li [l + l( v n)/ l 'n\ ■ Therefore, Eq. ( p2|) is still mean- 
ingful in the case of observables that do not depend on 
momentum. 

For a given potential V(q), it is convenient to devise a 
characteristic energy scale e (e.g., the barrier height for 
a double-well potential, the well depth for physical po- 
tentials that vanish at infinity, etc.) and a length scale a 
(such that variations of V comparable to e occur on this 
length scale) and write V(q) — ev(q/a). In this way one 
better deals with the dimensionless coordinate x = q/cr. 
If x m is the absolute minimum of v(x), the harmonic ap- 
proximation (HA) of the system is characterized by the 
frequency ujq = -\/ew"(a; m )/mCT 2 ; a dimensionless cou- 
pling parameter g for the system can be defined as the 
ratio between the HA quantum energy level splitting Hujq 
and the overall energy scale e, 
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hujo 



'n 2 v"{x n 



(25) 



The case of weak (strong) quantum effects occurs when 
g is small (large) compared to 1. In the following appli- 
cation we shall make use of the dimensionless variables 
only, i.e., energies are given in units of e, lengths in units 
of a, frequencies in units of u>q, and so on; the reduced 
temperature is t = l/(e/3). 

Let us then consider the double-well quartic potential 



2\2 



v(x) = (l~x 2 ) 



(26) 



v c g(x) = (1 — x ) — 3a (x) + 1 In n(x) 



It has two degenerate symmetric minima in x m = ±1, 
with v"(x m ) = 8. From the PQSCHA equations (j2|) 
and the definitions (^3|) and (|l6| ) we obtain 

(27) 
(28) 

(29) 
(30) 

(31) 
(32) 



v{x) = n 
/ 2 w 



n=l 

2 00 



(7m) 2 + f 2 (x) + 7rn7(n) 
(7m) 2 

[3x 2 + 3a{x) - 1] , 



a(x) 



_£_ 
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where 



7(n) 



y t 

t ^ (7m) 2 + f 2 (x) + 7rn7(n) 



j3hu)(x) g u>(x) 
2 2t ujq 



2t 



u 



and v n — (2-Kt/g)ujQn . Equations ( P9| ) and (pOj) have to 
be solved self-consistently. This task is done numerically; 
exact reference data can be obtained only for 7 = by 
numerical solution of the Schrodinger equation. 

In Figs. |] and || we report the shapes of the coordinate 
probability distribution V(x) = (S(x — x)) in the case of 
Ohmic damping, j(fn) = Two = const, at a very strong 
value of the coupling, g = 5; this gives a (nondissipative) 
ground state energy eo = 1.394, and the first excited level 
is ei = eo + 2.355. When dissipation is switched on V(x) 
tends towards the classical distribution V c ~ e - "^'*. 

Figure [| shows typical results found for the average po- 
tential energy v(t) = (v(x)). By comparing with the ex- 
act data at r = it appears that the PQSCHA gives very 
accurate results, in spite of the strong coupling. At lowest 
temperatures the PQSCHA tends to the ordinary self- 
consistent harmonic (or one-loop) approximation, since 
the effective Boltzmann factor tends to a S function. 

However, a more physical model should involve a non- 
Markovian memory damping function; in such a model 
there would be at least one more characteristic frequency 
scale (e.g., the frequency uju in the Drude model jjj) 
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above which 7(2) rapidly vanishes. In this case all aver- 
ages make sense, and whereas the averages of coordinate- 
dependent observables tend again to the classical behav- 
ior (in other words, the environment quenches the pure- 
quantum coordinate fluctuations), those of momentum- 
dependent ones go in the opposite direction due to the 
momentum exchanges with the environment. Of course, 
the PQSCHA expressions are even more useful in these 
physical situations. Further details and applications, as 
well as the extension to the case of many degrees of free- 
dom, will be given in a forthcoming paper. 

Useful discussions with Professor Ulrich Weiss are ac- 
knowledged. We are also grateful to Dr. Klaus Kirsten 
(University of Leipzig) for providing us with an elegant 
derivation of the tricky expansion of Eq. pa) . 




FIG. 1. Configuration density V(x) — {S(x — x)) of the 
double-well quartic potential for g = 5, t = 1, and different 
values of the Ohmic damping parameter r = 7 /loq . The filled 
circles are the exact result for F — 0; the dotted curve at 
r = 00 corresponds to the classical limit. 
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FIG. 3. Average potential energy v = (v(x)) of the dou- 
ble-well quartic potential vs. temperature, for g — 5 and 
different values of T = j/ujq. The filled circles are the exact 
result for r = 0; the dotted curve at V = 00 corresponds to 
the classical limit. 
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FIG. 2. Same as figure 1, for g = 5, t = 2. 
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